\[G = \prod_{i = 1}^k g_i\]

\[g_i = \frac{|4X_i - 2| + a_i}{1 + a_i}\]

Sensitivity Indices: analytical results

\[V_i = \frac{1}{3(1 + a_i)^2}; \]

\[ S_i = \frac{V_i}{\prod_{i=1}^k \left( V_i + 1 \right) - 1} \]

\[S_{Ti} = V_j\cdot\prod_{i\neq j}^k(1 + V_i)\]

Analytical sensitivity indices calculation

The function for calculating the variances and hence the sensitivity indices are defined below:

options(scipen = 999)

set.seed(100)

# Define the function to calculate V_i
calculate_Vi <- function(a_i) {
  return(1 / (3 * (1 + a_i)^2))
}

# Define the function to calculate S_i
calculate_Si <- function(Vi, V_product) {
  return(Vi / (V_product - 1))
}

# Define the function to calculate V_Ti
calculate_VTi <- function(Vi, V) {
  return(Vi * V)
}

# Define the function to calculate S_Ti
calculate_STi <- function(VTi, V_product) {
  return(VTi / (V_product - 1))
}

# Main script
a <- c(0, 2, 10, 30, 99, 99) # Example values for a
k <- length(a)

# Calculate V_i for each a_i
Vi <- sapply(a, calculate_Vi)

# Calculate the product of (1 + V_i) for total variance
V_product <- prod(1 + Vi)

# Calculate V_Ti for each a_i
VTi <- sapply(1:k, function(i) calculate_VTi(Vi[i], prod(1 + Vi[-i])))

# Calculate first-order sensitivity indices S_i
Si <- Si_analytical <- sapply(Vi, calculate_Si, V_product = V_product)

# Calculate total-order sensitivity indices S_Ti
STi <- STi_analytical <- sapply(1:k, function(i) calculate_STi(VTi[i], V_product))

The arbitrary coefficients are selected. The higher the value, the less significant the input parameter’s impact on the output.

a <- c(0, 2, 10, 30, 99, 99) # Example values for a
# a <- c(6.42,6.42,6.42,6.42,6.42,6.42)
k <- length(a)

The results are then printed:

# Output the results
cat("First-order sensitivity indices (Si):\n")
print(round(Si, 5))

cat("\nTotal-order sensitivity indices (STi):\n")
print(round(STi, 5))

The sum must be equal to 1.

Numerically calculate indices using the sensobol package.

  1. Define the sensitivity parameters and the function under inquiry:
set.seed(2)

k <- length(a)
params <- paste("x", 1:k, sep = "")
R <- 10^3
type <- "norm"
conf <- 0.95

sobol_Fun <- function(X, vector) {
  y <- 1
  
  for (j in 1:length(a)) {
    y <- y * (abs(4 * X[, j] - 2) + vector[j])/(1 + vector[j])
  }
  return(y)
}
  1. Sensitivity indices from sensobol are computed for each increasing power of 2 with lapply function:
N_values <- 2^(7:15)

results <- lapply(N_values, function(N) {
  mat <- sobol_matrices(N = N, params = params, type = "R")
  y <- sobol_Fun(mat, a)
  ind <- sobol_indices(Y = y, N = N, params = params, boot = TRUE, type = "norm", conf = 0.95, R = R)
  
  list(N = N, 
       Si = ind$results %>% filter(sensitivity == "Si") %>% pull(original), 
       STi = ind$results %>% filter(sensitivity == "Ti") %>% pull(original))
})
  1. the results object is converted to two usable data.frame objects, for an easier use:

The convergence plots for each input sensitivity can be plotted:

# Transform the data to be properly plotted
plot_data <- final_results_Si %>%
  pivot_longer(cols = -c(coefficients, Si_analytical),
               names_to = "N",
               values_to = "Si_value",
               names_prefix = "Si_") %>%
  mutate(N = as.numeric(N))

# Create the plot
Si_ggplot <- ggplot(plot_data, aes(x = N, y = Si_value, color = factor(coefficients))) +
                    geom_line() +
                    geom_point() +
                    geom_hline(aes(yintercept = Si_analytical, color = factor(coefficients)), 
                               linetype = "dashed") +
                    scale_x_log10(breaks = unique(plot_data$N)) +
                    labs(title = "Convergence of Si values to analytical values (dashed lines)",
                         x = "Sample size (N)",
                         y = "",
                         color = "Coefficients:") +
                    theme_minimal() +
                    theme(legend.position = "right")
                  
ggplotly(Si_ggplot)
# Transform the data to be properly plotted
plot_data <- final_results_STi %>%
  pivot_longer(cols = -c(coefficients, STi_analytical),
               names_to = "N",
               values_to = "STi_value",
               names_prefix = "STi_") %>%
  mutate(N = as.numeric(N))

# Create the plot
STi_ggplot <- ggplot(plot_data, aes(x = N, y = STi_value, color = factor(coefficients))) +
                     geom_line() +
                     geom_point() +
                     geom_hline(aes(yintercept = STi_analytical, color = factor(coefficients)), 
                                linetype = "dashed") +
                     scale_x_log10(breaks = unique(plot_data$N)) +
                     labs(title = "Convergence of STi values to analytical values (dashed lines)",
                          x = "Sample size (N)",
                          y = "",
                          color = "Coefficients:") +
                     theme_minimal() +
                     theme(legend.position = "right")

ggplotly(STi_ggplot)

Pass the cursor over the dots the read the sensitivity values.

Given the simplicity of the function, the convergence is reached quickly, even with a small sample size. These plots can thus be readily used for more complex functions, where convergence might be slower and more interesting to assess.

LS0tDQp0aXRsZTogIlNvYm9sIEcgd2l0aCBrID0gNiINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQoNCmBgYHtyIGluY2x1ZGU9RkFMU0V9DQpsaWJyYXJ5KHNlbnNvYm9sKQ0KbGlicmFyeSh0aWR5dmVyc2UpDQpsaWJyYXJ5KGdncGxvdDIpDQpsaWJyYXJ5KHBsb3RseSkNCmBgYA0KDQokJEcgPSBccHJvZF97aSA9IDF9XmsgZ19pJCQNCg0KJCRnX2kgPSBcZnJhY3t8NFhfaSAtIDJ8ICsgYV9pfXsxICsgYV9pfSQkDQoNCiMjIFNlbnNpdGl2aXR5IEluZGljZXM6IGFuYWx5dGljYWwgcmVzdWx0cw0KDQokJFZfaSA9IFxmcmFjezF9ezMoMSArIGFfaSleMn07ICQkDQoNCiQkIFNfaSA9IFxmcmFje1ZfaX17XHByb2Rfe2k9MX1eayBcbGVmdCggVl9pICsgMSBccmlnaHQpIC0gMX0gJCQNCg0KJCRTX3tUaX0gPSBWX2pcY2RvdFxwcm9kX3tpXG5lcSBqfV5rKDEgKyBWX2kpJCQNCg0KIyMjIEFuYWx5dGljYWwgc2Vuc2l0aXZpdHkgaW5kaWNlcyBjYWxjdWxhdGlvbg0KDQpUaGUgZnVuY3Rpb24gZm9yIGNhbGN1bGF0aW5nIHRoZSB2YXJpYW5jZXMgYW5kIGhlbmNlIHRoZSBzZW5zaXRpdml0eSBpbmRpY2VzIGFyZSBkZWZpbmVkIGJlbG93Og0KYGBge3J9DQpvcHRpb25zKHNjaXBlbiA9IDk5OSkNCg0Kc2V0LnNlZWQoMTAwKQ0KDQojIERlZmluZSB0aGUgZnVuY3Rpb24gdG8gY2FsY3VsYXRlIFZfaQ0KY2FsY3VsYXRlX1ZpIDwtIGZ1bmN0aW9uKGFfaSkgew0KICByZXR1cm4oMSAvICgzICogKDEgKyBhX2kpXjIpKQ0KfQ0KDQojIERlZmluZSB0aGUgZnVuY3Rpb24gdG8gY2FsY3VsYXRlIFNfaQ0KY2FsY3VsYXRlX1NpIDwtIGZ1bmN0aW9uKFZpLCBWX3Byb2R1Y3QpIHsNCiAgcmV0dXJuKFZpIC8gKFZfcHJvZHVjdCAtIDEpKQ0KfQ0KDQojIERlZmluZSB0aGUgZnVuY3Rpb24gdG8gY2FsY3VsYXRlIFZfVGkNCmNhbGN1bGF0ZV9WVGkgPC0gZnVuY3Rpb24oVmksIFYpIHsNCiAgcmV0dXJuKFZpICogVikNCn0NCg0KIyBEZWZpbmUgdGhlIGZ1bmN0aW9uIHRvIGNhbGN1bGF0ZSBTX1RpDQpjYWxjdWxhdGVfU1RpIDwtIGZ1bmN0aW9uKFZUaSwgVl9wcm9kdWN0KSB7DQogIHJldHVybihWVGkgLyAoVl9wcm9kdWN0IC0gMSkpDQp9DQoNCiMgTWFpbiBzY3JpcHQNCmEgPC0gYygwLCAyLCAxMCwgMzAsIDk5LCA5OSkgIyBFeGFtcGxlIHZhbHVlcyBmb3IgYQ0KayA8LSBsZW5ndGgoYSkNCg0KIyBDYWxjdWxhdGUgVl9pIGZvciBlYWNoIGFfaQ0KVmkgPC0gc2FwcGx5KGEsIGNhbGN1bGF0ZV9WaSkNCg0KIyBDYWxjdWxhdGUgdGhlIHByb2R1Y3Qgb2YgKDEgKyBWX2kpIGZvciB0b3RhbCB2YXJpYW5jZQ0KVl9wcm9kdWN0IDwtIHByb2QoMSArIFZpKQ0KDQojIENhbGN1bGF0ZSBWX1RpIGZvciBlYWNoIGFfaQ0KVlRpIDwtIHNhcHBseSgxOmssIGZ1bmN0aW9uKGkpIGNhbGN1bGF0ZV9WVGkoVmlbaV0sIHByb2QoMSArIFZpWy1pXSkpKQ0KDQojIENhbGN1bGF0ZSBmaXJzdC1vcmRlciBzZW5zaXRpdml0eSBpbmRpY2VzIFNfaQ0KU2kgPC0gU2lfYW5hbHl0aWNhbCA8LSBzYXBwbHkoVmksIGNhbGN1bGF0ZV9TaSwgVl9wcm9kdWN0ID0gVl9wcm9kdWN0KQ0KDQojIENhbGN1bGF0ZSB0b3RhbC1vcmRlciBzZW5zaXRpdml0eSBpbmRpY2VzIFNfVGkNClNUaSA8LSBTVGlfYW5hbHl0aWNhbCA8LSBzYXBwbHkoMTprLCBmdW5jdGlvbihpKSBjYWxjdWxhdGVfU1RpKFZUaVtpXSwgVl9wcm9kdWN0KSkNCmBgYA0KDQpUaGUgYXJiaXRyYXJ5IGNvZWZmaWNpZW50cyBhcmUgc2VsZWN0ZWQuIA0KVGhlIGhpZ2hlciB0aGUgdmFsdWUsIHRoZSBsZXNzIHNpZ25pZmljYW50IHRoZSBpbnB1dCBwYXJhbWV0ZXIncyBpbXBhY3Qgb24gdGhlIG91dHB1dC4NCg0KYGBge3J9DQphIDwtIGMoMCwgMiwgMTAsIDMwLCA5OSwgOTkpICMgRXhhbXBsZSB2YWx1ZXMgZm9yIGENCiMgYSA8LSBjKDYuNDIsNi40Miw2LjQyLDYuNDIsNi40Miw2LjQyKQ0KayA8LSBsZW5ndGgoYSkNCmBgYA0KDQpUaGUgcmVzdWx0cyBhcmUgdGhlbiBwcmludGVkOg0KYGBge3J9DQojIE91dHB1dCB0aGUgcmVzdWx0cw0KY2F0KCJGaXJzdC1vcmRlciBzZW5zaXRpdml0eSBpbmRpY2VzIChTaSk6XG4iKQ0KcHJpbnQocm91bmQoU2ksIDUpKQ0KDQpjYXQoIlxuVG90YWwtb3JkZXIgc2Vuc2l0aXZpdHkgaW5kaWNlcyAoU1RpKTpcbiIpDQpwcmludChyb3VuZChTVGksIDUpKQ0KYGBgDQpUaGUgc3VtIG11c3QgYmUgZXF1YWwgdG8gMS4NCmBgYHtyfQ0KDQpgYGANCg0KIyMjIE51bWVyaWNhbGx5IGNhbGN1bGF0ZSBpbmRpY2VzIHVzaW5nIHRoZSBgc2Vuc29ib2xgIHBhY2thZ2UuDQoNCjEuIERlZmluZSB0aGUgc2Vuc2l0aXZpdHkgcGFyYW1ldGVycyBhbmQgdGhlIGZ1bmN0aW9uIHVuZGVyIGlucXVpcnk6DQpgYGB7cn0NCnNldC5zZWVkKDIpDQoNCmsgPC0gbGVuZ3RoKGEpDQpwYXJhbXMgPC0gcGFzdGUoIngiLCAxOmssIHNlcCA9ICIiKQ0KUiA8LSAxMF4zDQp0eXBlIDwtICJub3JtIg0KY29uZiA8LSAwLjk1DQoNCnNvYm9sX0Z1biA8LSBmdW5jdGlvbihYLCB2ZWN0b3IpIHsNCiAgeSA8LSAxDQogIA0KICBmb3IgKGogaW4gMTpsZW5ndGgoYSkpIHsNCiAgICB5IDwtIHkgKiAoYWJzKDQgKiBYWywgal0gLSAyKSArIHZlY3RvcltqXSkvKDEgKyB2ZWN0b3Jbal0pDQogIH0NCiAgcmV0dXJuKHkpDQp9DQpgYGANCg0KMi4gU2Vuc2l0aXZpdHkgaW5kaWNlcyBmcm9tIGBzZW5zb2JvbGAgYXJlIGNvbXB1dGVkIGZvciBlYWNoIGluY3JlYXNpbmcgcG93ZXIgb2YgMiB3aXRoIGBsYXBwbHlgIGZ1bmN0aW9uOg0KYGBge3J9DQpOX3ZhbHVlcyA8LSAyXig3OjE1KQ0KDQpyZXN1bHRzIDwtIGxhcHBseShOX3ZhbHVlcywgZnVuY3Rpb24oTikgew0KICBtYXQgPC0gc29ib2xfbWF0cmljZXMoTiA9IE4sIHBhcmFtcyA9IHBhcmFtcywgdHlwZSA9ICJSIikNCiAgeSA8LSBzb2JvbF9GdW4obWF0LCBhKQ0KICBpbmQgPC0gc29ib2xfaW5kaWNlcyhZID0geSwgTiA9IE4sIHBhcmFtcyA9IHBhcmFtcywgYm9vdCA9IFRSVUUsIHR5cGUgPSAibm9ybSIsIGNvbmYgPSAwLjk1LCBSID0gUikNCiAgDQogIGxpc3QoTiA9IE4sIA0KICAgICAgIFNpID0gaW5kJHJlc3VsdHMgJT4lIGZpbHRlcihzZW5zaXRpdml0eSA9PSAiU2kiKSAlPiUgcHVsbChvcmlnaW5hbCksIA0KICAgICAgIFNUaSA9IGluZCRyZXN1bHRzICU+JSBmaWx0ZXIoc2Vuc2l0aXZpdHkgPT0gIlRpIikgJT4lIHB1bGwob3JpZ2luYWwpKQ0KfSkNCmBgYA0KDQozLiB0aGUgYHJlc3VsdHNgIG9iamVjdCBpcyBjb252ZXJ0ZWQgdG8gdHdvIHVzYWJsZSBkYXRhLmZyYW1lIG9iamVjdHMsIGZvciBhbiBlYXNpZXIgdXNlOg0KYGBge3J9DQojIENvbnZlcnQgcmVzdWx0cyB0byBhIG1vcmUgdXNhYmxlIGZvcm1hdA0KZmluYWxfcmVzdWx0cyA8LSBkYXRhLmZyYW1lKA0KICBjb2VmZmljaWVudHMgPSBhLA0KICBTaV9hbmFseXRpY2FsID0gU2ksDQogIFNUaV9hbmFseXRpY2FsID0gU1RpDQopDQoNCmZvciAoaSBpbiBzZXFfYWxvbmcocmVzdWx0cykpIHsNCiAgZmluYWxfcmVzdWx0c1twYXN0ZTAoIlNpXyIsIHJlc3VsdHNbW2ldXSROKV0gPC0gcmVzdWx0c1tbaV1dJFNpDQogIGZpbmFsX3Jlc3VsdHNbcGFzdGUwKCJTVGlfIiwgcmVzdWx0c1tbaV1dJE4pXSA8LSByZXN1bHRzW1tpXV0kU1RpDQp9DQoNCmZpbmFsX3Jlc3VsdHNfU2kgPC0gZmluYWxfcmVzdWx0cyAlPiUgc2VsZWN0KGNvZWZmaWNpZW50cywgc3RhcnRzX3dpdGgoIlNpIikpICU+JSByb3VuZCg1KQ0KZmluYWxfcmVzdWx0c19TVGkgPC0gZmluYWxfcmVzdWx0cyAlPiUgc2VsZWN0KGNvZWZmaWNpZW50cywgc3RhcnRzX3dpdGgoIlNUaSIpKSAlPiUgcm91bmQoNSkNCmBgYA0KDQpUaGUgY29udmVyZ2VuY2UgcGxvdHMgZm9yIGVhY2ggaW5wdXQgc2Vuc2l0aXZpdHkgY2FuIGJlIHBsb3R0ZWQ6DQpgYGB7cn0NCiMgVHJhbnNmb3JtIHRoZSBkYXRhIHRvIGJlIHByb3Blcmx5IHBsb3R0ZWQNCnBsb3RfZGF0YSA8LSBmaW5hbF9yZXN1bHRzX1NpICU+JQ0KICBwaXZvdF9sb25nZXIoY29scyA9IC1jKGNvZWZmaWNpZW50cywgU2lfYW5hbHl0aWNhbCksDQogICAgICAgICAgICAgICBuYW1lc190byA9ICJOIiwNCiAgICAgICAgICAgICAgIHZhbHVlc190byA9ICJTaV92YWx1ZSIsDQogICAgICAgICAgICAgICBuYW1lc19wcmVmaXggPSAiU2lfIikgJT4lDQogIG11dGF0ZShOID0gYXMubnVtZXJpYyhOKSkNCg0KIyBDcmVhdGUgdGhlIHBsb3QNClNpX2dncGxvdCA8LSBnZ3Bsb3QocGxvdF9kYXRhLCBhZXMoeCA9IE4sIHkgPSBTaV92YWx1ZSwgY29sb3IgPSBmYWN0b3IoY29lZmZpY2llbnRzKSkpICsNCiAgICAgICAgICAgICAgICAgICAgZ2VvbV9saW5lKCkgKw0KICAgICAgICAgICAgICAgICAgICBnZW9tX3BvaW50KCkgKw0KICAgICAgICAgICAgICAgICAgICBnZW9tX2hsaW5lKGFlcyh5aW50ZXJjZXB0ID0gU2lfYW5hbHl0aWNhbCwgY29sb3IgPSBmYWN0b3IoY29lZmZpY2llbnRzKSksIA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGxpbmV0eXBlID0gImRhc2hlZCIpICsNCiAgICAgICAgICAgICAgICAgICAgc2NhbGVfeF9sb2cxMChicmVha3MgPSB1bmlxdWUocGxvdF9kYXRhJE4pKSArDQogICAgICAgICAgICAgICAgICAgIGxhYnModGl0bGUgPSAiQ29udmVyZ2VuY2Ugb2YgU2kgdmFsdWVzIHRvIGFuYWx5dGljYWwgdmFsdWVzIChkYXNoZWQgbGluZXMpIiwNCiAgICAgICAgICAgICAgICAgICAgICAgICB4ID0gIlNhbXBsZSBzaXplIChOKSIsDQogICAgICAgICAgICAgICAgICAgICAgICAgeSA9ICIiLA0KICAgICAgICAgICAgICAgICAgICAgICAgIGNvbG9yID0gIkNvZWZmaWNpZW50czoiKSArDQogICAgICAgICAgICAgICAgICAgIHRoZW1lX21pbmltYWwoKSArDQogICAgICAgICAgICAgICAgICAgIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJyaWdodCIpDQogICAgICAgICAgICAgICAgICANCmdncGxvdGx5KFNpX2dncGxvdCkNCmBgYA0KDQpgYGB7cn0NCiMgVHJhbnNmb3JtIHRoZSBkYXRhIHRvIGJlIHByb3Blcmx5IHBsb3R0ZWQNCnBsb3RfZGF0YSA8LSBmaW5hbF9yZXN1bHRzX1NUaSAlPiUNCiAgcGl2b3RfbG9uZ2VyKGNvbHMgPSAtYyhjb2VmZmljaWVudHMsIFNUaV9hbmFseXRpY2FsKSwNCiAgICAgICAgICAgICAgIG5hbWVzX3RvID0gIk4iLA0KICAgICAgICAgICAgICAgdmFsdWVzX3RvID0gIlNUaV92YWx1ZSIsDQogICAgICAgICAgICAgICBuYW1lc19wcmVmaXggPSAiU1RpXyIpICU+JQ0KICBtdXRhdGUoTiA9IGFzLm51bWVyaWMoTikpDQoNCiMgQ3JlYXRlIHRoZSBwbG90DQpTVGlfZ2dwbG90IDwtIGdncGxvdChwbG90X2RhdGEsIGFlcyh4ID0gTiwgeSA9IFNUaV92YWx1ZSwgY29sb3IgPSBmYWN0b3IoY29lZmZpY2llbnRzKSkpICsNCiAgICAgICAgICAgICAgICAgICAgIGdlb21fbGluZSgpICsNCiAgICAgICAgICAgICAgICAgICAgIGdlb21fcG9pbnQoKSArDQogICAgICAgICAgICAgICAgICAgICBnZW9tX2hsaW5lKGFlcyh5aW50ZXJjZXB0ID0gU1RpX2FuYWx5dGljYWwsIGNvbG9yID0gZmFjdG9yKGNvZWZmaWNpZW50cykpLCANCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbGluZXR5cGUgPSAiZGFzaGVkIikgKw0KICAgICAgICAgICAgICAgICAgICAgc2NhbGVfeF9sb2cxMChicmVha3MgPSB1bmlxdWUocGxvdF9kYXRhJE4pKSArDQogICAgICAgICAgICAgICAgICAgICBsYWJzKHRpdGxlID0gIkNvbnZlcmdlbmNlIG9mIFNUaSB2YWx1ZXMgdG8gYW5hbHl0aWNhbCB2YWx1ZXMgKGRhc2hlZCBsaW5lcykiLA0KICAgICAgICAgICAgICAgICAgICAgICAgICB4ID0gIlNhbXBsZSBzaXplIChOKSIsDQogICAgICAgICAgICAgICAgICAgICAgICAgIHkgPSAiIiwNCiAgICAgICAgICAgICAgICAgICAgICAgICAgY29sb3IgPSAiQ29lZmZpY2llbnRzOiIpICsNCiAgICAgICAgICAgICAgICAgICAgIHRoZW1lX21pbmltYWwoKSArDQogICAgICAgICAgICAgICAgICAgICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAicmlnaHQiKQ0KDQpnZ3Bsb3RseShTVGlfZ2dwbG90KQ0KYGBgDQoNClBhc3MgIHRoZSBjdXJzb3Igb3ZlciB0aGUgZG90cyB0aGUgcmVhZCB0aGUgc2Vuc2l0aXZpdHkgdmFsdWVzLg0KDQoNCkdpdmVuIHRoZSBzaW1wbGljaXR5IG9mIHRoZSBmdW5jdGlvbiwgdGhlIGNvbnZlcmdlbmNlIGlzIHJlYWNoZWQgcXVpY2tseSwgZXZlbiB3aXRoIGEgc21hbGwgc2FtcGxlIHNpemUuIFRoZXNlIHBsb3RzIGNhbiB0aHVzIGJlIHJlYWRpbHkgdXNlZCBmb3IgbW9yZSBjb21wbGV4IGZ1bmN0aW9ucywgd2hlcmUgY29udmVyZ2VuY2UgbWlnaHQgYmUgc2xvd2VyIGFuZCBtb3JlIGludGVyZXN0aW5nIHRvIGFzc2Vzcy4gDQo=